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I""! Abstract 

A long wave multi-dimensional approximation of shallow water 
' waves is the bi-directional Benney-Luke equation. It yields the well- 

^ known Kadomtsev-Petviashvili equation in a quasi one-directional limit. 

r]^ A direct perturbation method is developed; it uses the underlying con- 

^ servation laws to determine the slow evolution of parameters of two 

O space dimensional, non-decaying web-type solutions to the Benney- 

' ^ Luke equation. New numerical simulations, based on windowing meth- 

ods which are effective for non-decaying data, are presented. These 
simulations support the analytical results and elucidate the relation- 
ship between the Kadomtsev-Petviashvilli and the Benney-Luke equa- 
tions and are also used to obtain amplitude information regarding 
particular web solutions. Additional dissipative perturbations to the 
Benney-Luke equation are also studied. 
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(3 1 Introduction 

The modeling of the propagation of waves on the surface of water is a classi- 
CN cal problem and remains an active and important area of research. However, 

the mathematical formulation of the problem, also known as the Euler water 
^ ^ wave equations, presents a number of difficulties such as strong nonlinearity 

^ and an unknown location of the boundary that makes finding solutions chal- 

^ lenging. To deal with this problem, researchers have employed a number of 

approximate models to the full system of equations describing surface flow. 

In weakly nonlinear water waves whose depth is shallow with respect 
to the wavelength of the wave, a maximally-balanced multi-dimensional ap- 
proximation to the Euler water wave equations is the Kadomtsev-Petviashvilli 
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(KP) equation. There are two signs associated with the KP equation de- 
pending on surface tension (ST); these are termed the KP-I (large ST) and 
KP-II equations (smaU ST). Here we assume surface tension is sufficiently 
small such that the KP-II equation is the relevant approximation to the 
Euler water wave equations, cf. [6j. Further, in order to derive the KP 
equation, one also assumes that the height of the wave is small compared to 
the depth, and that the flow is quasi one-dimensional; i.e. transverse varia- 
tions are slower than those in the primary direction of propagation (c/. [1] 
for a derivation and details). While other two-dimensional, shallow- water 
approximations to the Euler equations can be found, the KP equation is the 
unique model that is maximally-balanced. This, and the fact that it has a 
Lax Pair and many known closed form solutions makes the KP equation an 
important model. 

In recent years, a novel class of solutions, called 'web-solutions', of the 
KP equation have been discovered and analyzed. Building on work in |19j 
and the resonant solutions found in cf. |15] . |16j . in [10] the web solutions 
were obtained in terms of Wronskians via Sato's formalism |17] . Important 
properties and pertinent theorems concerning web-solutions can be found in 
Likewise, experimental and observational evidence exists, cf. [T], [14j, 
[22], that indicates that these web-solutions are present and persistent in 
nature. 

From a mathematical perspective it is important to understand whether 
these solutions are robust to perturbations. In this regard, it is useful to 
analyze improved approximations to the Euler equations which themselves 
yield the KP equation. An approximation that is the object of study in this 
paper is the Benney-Luke (BL) equation, cf. [8j. The BL equation is also a 
shallow water approximation to the Euler equation which arises naturally as 
an approximation to water waves. However, unlike the KP equation, the BL 
equation allows for two-directional waves. Also important in our study is 
an asymptotically equivalent modification of the BL equation, which unlike 
the KP equation, has a dispersion relationship in which the group velocity 
of arbitrarily large wave numbers is bounded; we consider this a BL type 
equation as well. The latter BL equation, aside from being a two-directional 
approximation to the Euler equations, is expected to be well posed over a 
broader class of initial conditions than the KP equation. 

Previous studies of the Benney-Luke equation have appeared in the liter- 
ature. The analysis in |15j begins with the BL equation, and a variant of the 
Benney-Luke equation was studied numerically in the context of the Mach- 
reflection problem in [13j . In each of these papers, particular web-solutions, 
so-called 'Y-junction' solutions, were studied. While this paper addresses 
the case of small surface tension in the BL equation, we note that the large 
surface tension case has also been studied. Rigorous results establishing 
the existence of lump soliton solutions to the BL equation have appeared, 
cf. [l8j, and the dynamics of these lump solitons have been investigated. 
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cf. M- 

Recently, we developed a perturbation theory of web solutions of the BL 
equation based techniques from the theory of integrable systems. We estab- 
lished in [3] that the web-solutions of the KP equation persist for asymp- 
totically long periods of time in the BL equation. This implies that the 
web-solutions are robust to the higher-order effects introduced by the BL 
equation. 

The purpose of the present paper is to: i) develop a direct perturbation 
method based on conservation laws of the BL equation and to ii) present new 
computations of the BL equation which are consistent with web type initial 
data initial data that is non-decaying. In one dimension the conservation 
law approach is widely used cf. j^. This is due to its inherent simplicity 
and does not require one to employ more complex integrable system meth- 
ods. However, in two dimensions since conserved quantities are expressed in 
terms of integrals, and the web-solutions do not decay in one of the spatial 
directions the integrals are infinite. To deal with this difficulty, we present 
a modified definition of the conserved quantities of KP that works for the 
web-solutions. 

With this modification, the conservation law approach now provides a 
direct method to perform perturbation analysis (via multiple scales meth- 
ods) on the web-solutions. As indicated above it does not require integrable 
systems machinery and readily achieves the same results. Using the con- 
venience of this conservation law approach, we are also able to analyze the 
impact of a typical type of dissipation on web-solutions in the Benney-Luke 
equation. We chose a linear local dissipative term which creates 'shelves' 
behind the web solution. In one dimension, due to mean term interactions, 
small amplitude shelves of long extent are known to arise from KdV models 
with such a dissipative term cf. |5J. 

While the type of dissipation we study is particular, the method we 
present can, in principle, be used for many other, more physically realistic, 
dissipation models, cf. [12] . In this paper our aim is not to decide on the 
best dissipative model. Instead, we explain how the direct conservation 
law approach needs to be modified when dissipation is present. We further 
exhibit the growth of shelves from the linear local dissipative model. This 
method can also be extended to higher order approximations of water waves. 
However, carrying out this effort is outside the scope of this paper. 

The asymptotic methods in this paper rely on separating the web-solutions 
into what we define as near and far field components. Using our conserva- 
tion law based asymptotic method, the far field components are studied 
analytically. The near field component is not required in the leading order 
perturbation analysis. However, the near field is where the various parts of 
a web-solution interact. Therefore, we employ new numerical simulations 
of the BL equation to investigate the near field, or interaction region. In 
this paper, we in particular look at both Y-junctions and another class of 
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web-solutions, the X- waves (or 0-type solutions in the terminology of 

We also examine the Y-junction solution further by studying the amplifi- 
cation of the so called "Mach stem" , cf. [16j . By varying the small parameter 
representing the 'shallowness' of the water and the angle of opening of the 
'Y', we show how the relative amplification (see Section [s]) of the Mach stem 
in the BL equation varies. We find that the relative amplitude, or amplifi- 
cation ratio, of the stem increases as the water becomes shallower and this 
result holds uniformly over several different angles. As the water becomes 
shallower, or as the BL equation tends to the KP equation, the amplification 
ratios of the BL equation tend increase to those of the KP equation. For 
the KP Y-junction solution, the largest amplification ratio is known to be 
four For the BL equation we find that this ratio is reduced by 0(e). 
An interesting open question is: how much is this ratio reduced in the full 
water wave equations? 

Thus the numerical study yields results consistent with asymptotic ar- 
guments and is also useful for predictive purposes. We also numerically 
solve the BL equation with our chosen dissipation term and show that small 
shelves of long extent form in the wake of the Y-junction; this is also con- 
sistent with previous one dimensional results, cf. [5j. All of these results are 
new in the literature. 

In summary, in this paper we: 

• Develop a new, conservation law based perturbation method applicable 
to any web solution. This significantly simplifies previous, integrable 
systems based approaches, and thus allows us to model perturbations 
to the KP equation more readily. 

• Examine, for the first time, the role dissipation has on web-solutions. 
While we have chosen a particular dissipation model, our methods are, 
in principle, applicable to any dissipation model. 

• Present new numerical results showing how higher order shallow wa- 
ter effects influence the evolution of the interaction region of web- 
solutions. In the case of Y-junctions, we do this over a large number 
of parameter values and investigate relative amplification ratios; these 
results are consistent with KP theory. 

While in this paper we only study web-solutions in the context of shallow 
water flow, we believe that the combination of asymptotic and numerical 
methods presented here will be useful in other nonlinear problems involving 
non-decaying proflles. 
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1.1 The Benney— Luke Equation and Web-Solutions of the 
KP Equation 

The Benney-Luke (BL) equation is given in non-dimensional form by [1] 

qtt -Aq + aeK^q + e(9f|Vg|2 + qtKq) = 0, (1) 

where small dispersion, slow transverse variation and weak nonlinearity are 
all balanced. These small effects are denoted by e, |e| <C 1; also q is the 
velocity potential, V = (^dx, e^^'^dy) , A = d'^ + edy, and a = a — t^, where 
a is related to the surface tension. Since the BL equation is a long wave 
approximation, as an initial value problem, it can and does suffer from 
arbitrarily large growth rates. 

By replacing (?| by we shall use a a regularized version of the BL 
equation which is asymptotically equivalent to the above BL equation ([T]) 
and is written in the form 

qtt -Aq + aeAqtt + €{dt\Vq\'^ + qtAq) = 0. (2) 

We use the transformation ^ = x — t, r = et, so that 

= dx, dt = -d^ + edr- 

In this paper we only consider small surface tension, hence o" ~ 0, so that 
a ~ —1/3 and thus this regularized BL equation is linearly well-posed; it 

has the dispersion relationship o;^ = ^|l||fc|2 ' ^^^re k = {kx,\/eky). 

Introducing the scaling A = {—8a)^^^, P = l"^ = ^ ~ ^^"^ 
letting q = ^0(/3^, 7?/, Jr), ^ becomes 

Defining the operator (see j2]) 

d^^uiC,y,t) = - J u{z,y,t)dz - - u{z,y,t)dz, 

and letting u = (j)^ and e = eA the BL equation we consider is given by 

(— 4mt- + u^^^ + 6uu^)^ + Suyy — ed^F{u) = 0, (3) 
where F{u) is given by 

F{u) = 2d^^UrT — '2,u^^r — ^u^yy — 39^(9^"^%)^ 
—4uUr — 3ud7^Uyy — 2u^d7^Ur- 
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Hereafter we replace e by e (also note that for small surface tension A ~ 1.2). 
The definition of d^^ is chosen so that if /, (7 G -^^^(K), then /j^ 
is well defined, and 

/ fiOdT'giOd^ = - [ 9(097' fiOd^. 

Alternatively the perturbed BL equation is written as 

Ur = K(u) - '-F[u) (4) 

where K(u) = |tt^5^ + \uu^ + \d'^^Uyy . 

As seen from the rescaled BL equation ([3]), the KP equation gives the 
leading order behavior of the BL equation. We now introduce the web- 
solutions to KP which serve as the particular leading order behavior of 
interest. A web solution, say w(^,y,T) to the KP equation in Wronskian 
form is given by (cf. [11]; note for the general N-soliton solution in Hirota 
form see [19j) 

w(i,y,T) = 2dl\og(n(i,y,T)), (5) 

where y, r) = Wr(/i, • • • , /tv), where Wr denote the Wronskian of the 
functions fi {cf. pTT])- The functions fi have the particular form 

M 

f, = ^b,,e'^, (6) 
i=i 

with bij constant, < M, and Oj{^,y,T) = kj^ + kjy + kjT, where ki < 
was proved in |llj that as y — t- 00, there are N lines, and as 
y — )• —00, M — N lines, in the ^,y plane with slopes Cij = —(ki + kj)~^, 
j > i, such that 

..~l(fc,-fe,)W(t^±^). (7) 

where the phases 9ij are constant in this pure KP solution. 

We now introduce a perturbation scheme for studying the behavior of 
solutions y, r) to the BL equation of the form 

u(^, y, t) = w(^, y, r, T) + es{^, y,T) ^ , 

where T = er. The leading order term ^i; is a web-solution to the KP equa- 
tion in which we allow the terms kj to vary in the slow time T. Motivated 
by the form of the web solution, it is convenient to separate the solution 
w into far field (\y\ » 1) and near field (\y\ = 0(1)) components (see also 
[3]). The far- field coordinates along a given ray are denoted as Xij and Yij, 
where j > i, and Xij is determined by the differential equations 

d^Xij = 1, dyXij = ki -\- kj J d^Xij = kj -\- k^ -\- kikj^ (8) 
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and Yij is given by 

Yij = y- (kj + ki)C, (9) 
where in the large \y\ hmit derivatives with respect to Yij are considered 
asymptotically negligible. We also require that ki{T) + kj{T) = Cij, where 
Cij is independent of r and T; hence dxXij = 0. This requirement means 
there is no secular growth in drw for large y. 
We assume along \Yij\ oo that 

u ~ w{Xij,T) + es{Xij,T) + ■■■ . 

In this limit, the equation for w is, 



dx,, ( -4^^^^^^x,, + d'x^.w + 6wwx,, ] + S{kj + kifd],^.w = 



where kij = kj — ki. Integrating once in Xij, we get 

-kfjWXij + d^..w + 6wwxij = 
which has the soliton solution 

and which is consistent with the previous asymptotic arguments. The func- 
tions kij = kij (T) are functions of the slow time and 9^'^^ (T) is an arbitrary 
phase; these functions are determined later in the perturbation scheme. 

At the next order, noting for convenience that we drop the subscripts on 
Xij, the equation for s{X,t) is 

-Asr - kfjsx + d\s + Q{ws)x = dxF{X) + Awt, 
It is also noted that 

F(X)=2S^.(X)- (^^(^. + fc,)^ + 3M^^ai. (11) 



-\2^\^ + "^{kj + kif\w\X), 

fvij Z 



and 

Wt 



+ f ^(X + ^g)) + dre^) wx. (12) 

Kij \ Kij J 

Since dxF[X) + Awt is independent of the fast time r, we separate the 
perturbation s into s = h{X) + g{X,T), so that 

-kfjhx + d\h + Q{wh)x =dxF{X) + 4.wt, (13) 
-Agr - kfjgx + d^x9 + G{wg)x =0, 
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with g{X,0) = -h{X). 

The slowly varying phases 6^^^ {T) were determined in [3] via the addi- 
tional orthogonality condition 

f w{X)h{X)dX = 0. (14) 
Jr 

which removed unbounded growth as r — )• oo. We note that orthogonal- 
ity conditions such as these usually correspond to physical constraints, i.e. 
conservation laws. 

In Section 2, from the conservation law approach a direct and physically 
motivated derivation of the condition ( 14 ) is given; also in Section 2 the role 



of dissipation is examined. In Section 3, we present our numerical findings 
concerning the evolution of several different types of web-solutions. These 
results indicate that web-solutions are robust to these perturbative effects 
over asymptotically significant lengths of time. The numerical study of the 
amplification of the stem of a Y-Junction in the BL equation, and the effect 
of dissipation, is also presented in Section 3. The Appendix collects some of 
the technical details of the approach used in Sections 2 and 3. 

2 Conservation Law Scheme for Arbitrary Web- 
Solutions 

In order to determine the evolution of the slow coefficients, we use a method 
(see [5]), which employs the conserved quantities of the leading order inte- 
grable problem to determine the evolution of slow variables. We use conser- 
vation of energy which we define by the following iterated integral 

= U f u^dA = lim ^ r [ u\i,y,T)didy. 

The average along the y-axis is taken because the line solutions to the KP 
equation do not decay along the y-axis. It is also important to point out 
that for web-solutions and their perturbations that 

lim / y / u^{ty,r)dyd^ ^ [ lim y / u^{(,,y,T)dyd^, 

and so the placement of the limit is important. 

Using the above definition of energy it follows that 

dr£{u) = ( [uUr dn = ( I uK{u)dn " 7 ( / uF{u)d^ 

\Jr / y \Jm. I y 4 XJjg I y 

The analytic issues concerned with differentiating the energy functional are 
discussed in the Appendix of this paper. From the definition of -ftr(ti), we 
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expect that 



which would then give 



[ uK{u)dn = 0, 



(15) 



UUrd^ 



uF{u)di 



(16) 



Note, due to the average, in order to show (15), we must estabhsh that 



/ ud^ ^UyydS, ) = 0. 

JR / y 



(17) 



To show this, some care must be used to accommodate for the presence 
of the average which makes integration by parts more compHcated. This 
is worked out in the Appendix for web solutions to the KP equation. We 
further assume it holds for the perturbation u = w + es + ■ ■ ■ , so that taking 



( 15 ) to be true is valid to the order of the asymptotics. 

Since we are assuming u{^, y, r) is, in particular regions of the plane, a 
solution of the form 



u(e, y, r) = w{Xij,T) + es{Xij,T) + 0{e' 



(18) 



and is otherwise at most 0{e) outside the regions of interest, the ansatz for 
u implies Wr = cwt- Further, based on the one dimensional problem it is 
natural to assume that 



/ 

Jr 



WSrd^ 



oie), 



(19) 



(this is discussed in more detail below), then we get 

/ w{wT + ]F{w))dn = 0. 
\Jr 4 / 



(20) 



We now demonstrate, in some detail, how to use the solvability condi- 



tions associated with (13)) in order to reduce (|16| ) to (|19[ ) and (20). As 
|y| oo, using the coordinate change given by (IsFand ([9]), w depends only 
on Xij, hence becoming stationary with respect to r. The perturbation 
s{Xij,T) is written as s(Xjj,r) = h{Xij) + g{Xij,T) in region Rij (see Fig 
([l])) where h solves 



K'{w)h 



-F{w) + wt- 



(21) 



Note, (21) is the same as (13). Then in Rij, where we ignore the Yij coor- 
dinate and thus reduce the problem to one dimesnion, i.e. Xij, as shown in 
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Inleraetion 
Region / 



Figure 1: Far-Field Region 



the introduction, we have {K'{w)) w = 0. In order to ensure that a solution 



h exists to (21), the Fredholm alternative enforces to leading order that 

j w{wT + dXij = 0. (22) 

This result gives us immediately that the terms kij are independent of the 
slow time T. To show this, we first note that 



wiX,,)F{X,j)dXi. 



wiX,,)dx,,F{Xij)dX,j = 0, 



(23) 



since w is even in X, and from (11), so is F. Using (22) then shows that 



wwTdXij = 0. 



On the other hand, using (12), we have 



wwxdXi. 



3 drhj I „„2 
2 



w'^{Xij)dXij. 



Since J^w'^{Xij)dXij ^ 0, this implies that = 0. Having now shown 

that the coefficients kij do not vary slowly, and given our assumption that 
kj + ki does not depend on T, evidently the coefficients kj are also indepen- 
dent of the slow time T . 



We now present an argument that establishes (20) using (22). First, we 
have the identity 

^ / F{w)\ f ( F{w)\ 

+ ,^'7^»(„. + £M)« 

/ wiwT + ^^^j^didy. 
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Taking S> 1, and letting L — )• oo, we can write the integrals 



w[wT^ -7-^ dydi, I w [wT-\ dyd^ 



La 



as sums over the regions Rij since in between the regions, w and its deriva- 
tives are exponentially small. Then, over any region Rij, we have 



W Wt + 



F{w) 



d^dy ~ Jj. 



Ya 



w Wt + 



F{w) 



dXj jdYj 



where Jij = 1 + (fc; + kj)~'^ is the Jacobian of the coordinate transformation 
and Ya is some constant value of Yij that only depends on the choice of La- 
Using (22) then establishes that 



w [ wx H — ) dt,dy ~ 0. 



Given that 



Ri 



1 Z"^" /■ / Fiw) 
lim - / w{wT + ] didy = 0, 

L^oo L J_r Ju \ 4 



since La is constant, we finally have 



1 



L 



hm — / I w [wt + . 



F{w) 



d(,dy ~ 0, 



or that (20) holds 



We now turn to determining the phases 6f^\T). To do this, using (20) 



(16) becomes, noting also that s(^, y, 0) = 0, 



swd^ 



0. 



(24) 



This global relation can be separated into near and far-field components; 
note the integrals in the interaction region [—La, La] are negligible in com- 
parison with the far field. Then in the far-field, we can write this integral 
over regions Rij (again see Figure [l] for clarification). If we assume that 
after a short time the dominant part of s{Xij,T) is given by h{Xij), then 



since kj is constant wt = 9x0^-^ wx^, so that (|l3l) becomes 



-kfjh + d\^^h + 6wh = F{Xij) + drefj^w. 



One gets from the global condition (24) that in Rij, again ignoring Yij, 



h{X,j)w{Xij)dXij = 0. 



(25) 
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This is exactly the orthogonality condition derived in [H], i.e. (14), 
that was used to determine the evolution of the difference between phases. 
To compute this phase evolution, first introduce the transformation X = 
kij{Xij + 9^j^)/2, and then let q{X) = 2sech^(X), which turns (13) after one 
integration in X into 

-4/i + + 6g/i = cjjq + cf^q"^ + cfj5|g 

where the coefficients (^j,P = li2,3, which implicitly depend on the region 
of interest Rij, are given by 



This leads to the solution 

h{X) = -2q{X) ^(-l + Xtanh(X)) - ^ + ^Xtanh(X) . 

\ 8 2 / 



Using (g gives c\j = |(-2c2j./3 + cfj/2) or 

{kf - 28kfkj - 54k^k^ - 28kik^ + fcf) 
drOlf = — — ^ ^ ^. (26) 

With this, we now can compute the speed, Vij, of each part of the far field 
via the formula 

Vi,j = drXi,+edTefj, (27) 

where drXij is given in ([s]). 

Therefore, using a conservation law approach, we have seen how the 
results of [3] can be derived in a way that makes no use of integrable systems. 
This makes the perturbation method presented in this paper more broadly 
applicable than that presented in ^ . 



2.1 Dissipation 

The Benney-Luke equation conserves energy, and is a Hamiltonian system 
(see [H]). Of course, real physical systems such as water waves are not 
perfectly conservative. So it is useful and instructive to study the effect of 
adding typical damping terms to the conservative model. While the dis- 
sipative model we use is special, it is chosen to show how to deal with 
methodological difficulties created by dissipation models. It is also instruc- 
tive to study this dissipation model since it causes the formation of small 
amplitude shelves with long extent. 
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For the Benney-Luke equation, we add a linear local term which leads 
us to study an equation of the form 

—4ur + u^^^ + 6uu^ + ^d^^Uyy — eju — eF{u) = 0, 

where the Fj are dispersive terms and the constant 7 > represents the 
magnitude of damping. We can then define the function Fi(u) = 7U + F{u) 
and repeat the analysis from above. Matching terms of ©(e) after expanding 
dr£{u) gives 



/ wwTd^/ + dr ( f wsd^\ = — - 

Jr I y \Jr I y 4 



(7'u; + wF{w))d(^ 

' y 



and repeating the Predholm alternative argument that led to ( 20 ) , yields 



( / wwrdn =-\{ [ {-fw^ + wF{w))dn . (28) 

\Jr I y 4 I y 



Now we seperate (28) over the regions Rij. We showed in the previous 



section, i.e. Equation (|23|), that 

wF{w)dXij = 0. 



Then, if we naively use the asymptotic form of the solution given by ( 10 ) 
and the corresponding representation for wt in we get 



2 kij 



[ w\Xi,)dX,j = -7 / w\Xij)dXij, 



so that we have 



dkjj _ _7 
dT ~ 6 



Coupling this condition with the requirement that ki{T) + kj{T) = dj 
creates a difficulty. For example, if we have ki, k2, and k^, then a web- 
solution whose slopes are determined by say ki and k2, and ki and k^, would 
have two different values for k^{T). Thus we see that the ansatz given by 



equation ( 10 ) for the leading order behavior is not appropriate globally. 

To deal with these inconsistencies, we use a modified ansatz for w (the 
index i, j on w and ■0;, I = 1,2,3, is understood) of the form 



w = sech 



where t]{T) = {kj{0) - A;i(0))e-'^^/^ and 

d^X = 1, dyX = kj{0) + k,{0), drX = V? + + i^l 
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with 



MT) = \ (MO) + m - rjiT)), MT) = I {k,{0) + h{0) + r^{T)). 



Note, the functions rj{T), 'ipi{T), and ip2{T) vary in each region Rij. 

To deter 
relationship 



To determine again noting that s(^, y, 0) = 0, we employ the global 



wsdn =0. 
' y 

To make use of this identity, we first assume that on Rij that the dominant 
contribution to s is given by the solution to the stationary equation 

-rj'^hx+d^^h + 6{wh)x = dxF{X) + iwx + -fw, (29) 



which is the same equation as (13) except now taking the dissipation into 
account. Using the transformations from the previous section we solve for 
h and derive the phase equation 

^^^^^ =6^ + 24 • 

3 Computation of the Benney-Luke Equation 

We first note that the Benney-Luke equation is second order in time. It 
is convenient for the numerics to transform to a system by introducing the 
variable v = Ur- In order to simulate the Benney-Luke equation numerically, 
we use a windowing method developed in |20j . This method introduces a 
smooth function of compact support, say V{y), such that V{y) is nearly one 
over some interval in y, say [—Ly + 5, Ly — 5], and V has support in [—Ly, Ly]. 
We use the function 

where eps is on the order of 10^^^ (approximately machine precision) and 
write the solution to the Benney-Luke equation as 

«).n.)(:) + a- >'(.))(: 

We use the leading order asymptotic solution computed in the previous 
sections to evaluate u and v in the far-field with distances greater than 
0{\Ly\) for times 0{l/e). We denote these far field solutions as Uasy and 

Vasy ■ 

We now wish to d-ctcriiiinc \^(^y^u. DcfinG th.G functions Unr — 
Vnr = Vv, SO that 

n ~ n„r + (1 - V)Uasy, f ~ t^nr + (1 - V)Vasy 
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Using the fact that Uasy satisfies the KP equation in the far field, and keeping 
only terms through 0(e), by substituting u and v into the BL equation, it 
then follows that Unr satisfies the equation 

where BC{u) denotes the Benney-Luke equation. The term Ff is given by 



2e V ^|(6(1 - V)UnrUasy - 3F(1 - V)ul,y) - QV' OyUasy - W'Uas, 

which has support only in the intervals [—Ly, —Ly+S\ and [Ly — 6, Ly\. Note, 
we ignore the window effects from the nonlinearities of 0{e) since these terms 
are smaller and isolated to have support in the region \—Ly, —Ly + 5] and 
[Ly — 5, Ly]. If we work on a large enough domain, the effect of the error 
should be nominal far from the boundary since the error propagates with 
finite speed. 

The functions Unr and Vnr satisfy at r = the boundary conditions 

UnriC, -Ly, 0) = Unr{^, Ly, 0) = 0. 

At later times, we enforce periodic boundary conditions in the y-coordinate 
and choose the window large enough to ensure that 

Unri^, -Ly, r) = UnriC, Ly, t) ~ 0. 

Given the exponential decay in ^, we also restrict the domain in ^ to the 
interval [—L^,L^] and enforce periodic boundary conditions in ^ and y in 
order to numerically solve for Unr- This gives us the advantage of being able 
to use pseudo-spectral methods. 

However, we have to find the operator d^^ in terms of Fourier series. 
Writing a typical periodic function, say u, as 

oo 
j,l=-oo 



we calculate the inverse derivative of u via 



^ m ^ ^ I m ^ ^ ' I 



where wc have used d^"^ = f^i^ — \ f^^ ■ In the Appendix, wc show that 
the pseudo-spectral representation of the (j, /) mode of the inverse derivative 
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is given by 

(d^'n)], =(1 - 6,o)^ a,,(r) - 5,o^ £ ^ 

k=-N+l,k=/=0 

where aji{T) = (— l)'+-^aj/(r) and Nt is the number of modes used in the 
pseudo-spectral approximation. 

3.1 Evolution of the Y- Junction and X-Wave in the Benney— 
Luke Equation 

The following figures show top-down surface plots of the numerical evolution 
of various web-solution profiles. The time stepping algorithm used was the 
ETDRK4 method [21], which is a variant of the 4th order Runge-Kutta 
method. 




(a) Long Stem, = 10"**, e = .2 (b) Long Stem, fcg = 10"^, e = .1 

Figure 3: Long Stem X- Waves: r = 1/e 

In each figure, the initial conditions were chosen as 'w{^, y, 0) and Wt{S,, y, 0) 
as defined earlier, see ^ and (|6j). The X -wave solutions in Figures [2]and|3] 
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>. 




50 ^0 50 

(a) Y Junction: k = 1, e = .2 

Figure 4: Y- Junctions: k = 1, r = l/e 



^0 50 

(b) Y Junction: fc = 1, e = .1 



correspond to the classical two-soliton solutions discussed in [6] and [16] . In 
the language used earlier, we choose the matrix of coefficients B with entries 
Oij as, see (§ and (§, 



B 



1 6i2 
1 624 



This gives the function y, r) as 

n =(^3 - ki)e^'+^'' + {k4 - ki)b24e^'^^^ 

+ ih - k2)h2e^^+^' + {ki - k2)bub24e^'^^^. 



We then choose bu = l^Zll , &24 = and we set ki = -k^, k2 = -k^, 

> 0, and k^ = 1 + ^3. Given this choice of coefficients bij and ki, by 
letting fcs — 7- 0, one can let the length of the stem grow as shown in Figure 
|3] Note that the radiation is more significant in Figures [2^, [3^ where e = 0.2 
than Figures [2|3, [3|3 where e = 0.1. 

A different case can be investigated by setting k^ = k2- Then Q becomes 



n = {k2- ki) [ e 
If we then choose 



6*1 +6*2 _^ ^4 - fci 
k2 - ki 

k2 - ki 



k2 - ki 
ki — ki 

ki — ki' /c4 — A;2 ' 

and relabel 64 as ^3, we then get Q = (/c2 — ki)Q, where 



^24 
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We note that 



a|ln((/c2 - A;i)f]) =a|ln(j^), 
and so from this point on we use w = 25| ln(0) as a solution to the KP 



equation. The function Q corresponds to the B matrix 

B = 



1 -1 
1 1 
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which shows that the X-wave has passed toaM = 3, N = 2 web-solution. 
This is an example of the Y-junction. We choose the parameters kj such 
that 

ki = —, k2 = -^—1 and = (32) 

where < A; < 1. Again, we note that the radiation is more significant in 
Figure where e = 0.2 than Figure [Ija where e = 0.1. 

These Y-junctions have a ray along y — t- — oo with maximum amplitude 
l^maxl = ^^~^2^ ; another ray along the angle arctan(— /c) of amplitude 1/2, 
and finally a ray of amplitude k'^/2 moving at an angle of 7r/4 where the 
angles are measured relative to the y-axis. Thus as A; — )• 1, the KP equation 
predicts a four-fold amplification of the ray along y — t- — oo relative to the 
amplitude of the ray along the angle arctan(— A;). Likewise, as we vary k, 
we move between a one-dimensional line soliton solution when k=0 up to a 
Y-junction with maximal amplification ratio of a factor four when k = 1, 
see Figure [4| Alternatively we can use the formalism in [19] to get these 
initial conditions. 

The figures represent numerical simulations run on time scales long 
enough to allow the next order terms that distinguish the BL equation from 
the KP equation to have a significant impact. This effect manifests itself in 
weak dispersive tails indicated in the figures. However, in all six figures, it 
can be seen that the primary wave remains essentially unchanged over time 
scales that are the reciprocal of the order of magnitude of the perturbation. 
We also point out that the figures show that as e decreases, or as the water 
becomes shallower, the amount of radiation decreases. These figures indi- 
cate that the web-solutions studied here are stable with respect to evolution 
in the BL equation. Hence this result provides quantitative evidence that 
the class of web-solutions examined in this paper are robust with respect to 
the perturbations introduced by the BL equation. 



3.2 Amplification Ratio for the Y- Junction in the Benney— 
Luke Equation 

The problem of the amplification ratio of the stem in a Y-junction is an 
interesting issue that was first studied in [16j in the context of the KP 
equation. We now discuss how the BL equation affects the amplification 
ratio, and how this class of KP solutions behaves under the infiuence of the 
BL equation. 

As indicated above, taking a Y-junction with parameter values given by 



(32), and its corresponding derivative with respect to time, as Cauchy data. 
Figure [s] shows the amplification ratio (vertical axis) of a Y-junction in the 
BL equation for different values of e and k. For each curve, the simulations 
were run for times up to r = ^, thus allowing time for the BL equation 
to have an asymptotically significant impact. The domain size was chosen 
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to be = 50 and Ly = 48 with 512 modes used for the pseudo-spectral 
approximation. The size of the rays of the Y-junction were measured at 
Ly = ±30 so as to avoid any effects from the windowing. In order to find 
the wave amphtude, we averaged the height of the wave over seven mesh 
points around the grid point corresponding to Ly = ±30; this minimizes the 
impact of any spurious oscillations due to the numerical method. Finally, 
we chose k = .1, • • • ,1 in intervals of a tenth. Figures 4a and 4b are plots 
of the case A; = 1 for e = .2 and e = .1 respectively. 




Figure 5: Amplification Ratio of the "Mach Stem" in the BL equation. 



We also note that the phase, or speed correction, along the ray y 
— oo is given by dxO^^ = ^^~^q^ so that the speed, see (27), of the ray in 
the far field should be — ( 1 ± e^^w^ )) where the negative sign indicates 



96 

the ray moves to the left. While this correction goes into the windowing 
approximation, as a consistency check, we also numerically computed the 
speed of the ray at y = —30 for e = .1 and k = 1. After averaging to 
take into account errors from discretization, we got a computed speed of 
— 1.0090, while the predicted speed from the asymptotics is —1.0167, so the 
agreement is very good. 

Figure [5] indicates that the numerically generated solutions tend to the 
exact results for the KP equation as e decreases and shows that increasing 
e decreases the amplification ratio for each k. As indicated above, we note 
that larger values of e introduce more dispersive radiation. This also implies 
lower pulse heights due to energy conservation which in turn suggests that 
there will be a somewhat smaller amplification ratio than four at A; = 1 of 
0(e) {i.e. about 3.9 for e = .1), but still much larger than linear theory. 
This result also agrees with the fact that the BL equation contains the KP 
equation as an asymptotic limit. 

The variation of the amplification as a function of e holds uniformly 
for all k. This indicates that the numerics is well-behaved as we vary the 
shape of the initial conditions. Therefore, we find that the combination of 
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asymptotic evaluation of the far-field, pseudo-spectral method, ETDRK4, 
and windowing is an effective way to approach problems with non-decaying 
solutions. 

3.3 Dissipation in the Y-Junction 

We next provide numerical results corresponding to the linear dissipation 
model discussed earlier. As a leading order solution, we use the Y-junction 
from the previous section with k = 1 so that /ci3 = 2, and we set the 
coefficient of dissipation 7 = 2. However, we also take into account the 
dissipation in the far- field using the asymptotic method shown earlier; this 
is how we find Uasym and Vasym- The size of the domain is = 50 and 
Ly = 48, with 512 modes used in the pseudo-spectral method. We measure 
the maximum height of the numerical solution in the region —30 < y < 30 
in order to avoid any effects from the windowing. The plot in Figure [6] gives 
a log plot of the amplitude decay. 

Thus, at r = 1/e, the asymptotic theory predicts that the maximum 
amplitude of the ray should be 2 • e~'^^ = 1.033, while the numerics gives 
the maximum amplitude as 2 • e~'^^ = 1.15. Thus the numerics confirms 
the asymptotic prediction since the discrepancy between the two results is 
0{e). We also look at a cross-section of the ray at y = —30 in Figure p| As 



1.5 




(a) Amplitude Decay Comparison (b) Solution profile at y = — 30 



Figure 6: 7 = 2, e = .1 

shown in the figure, to the right of the peak a shelf forms. This is consistent 
with the dissipative theory of KdV for linear dissipative models of wave 
amplitude |5J ; this is another novel feature of web-solutions not yet reported 
in the literature. 
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Appendix 

Basic Properties of the Average 

While this paper does not present rigorous results, we include some analyti- 
cal details that add support to the formal arguments presented in the paper. 
Supposing /(^, y) e L^{Rx [—L, L]), which is to say that /(^, y) is in L^(M) 
with respect to ^ and L^{[—L,L]) in y, then the quantity 

/ f{i,y)dydi 

-L Jr 

is well-defined, and by the Fubini-Tonelli theorem, and assuming the limits 
exist, we have that 

lim I J r /(^,y)«= lim y f{U)d^dy. 

We also have that if f{^,y) and are in L^(M x [—L,L]), then 

7 mr)d^) = lim y T / M^,y)d^dy = 0. 

\J'R I y '-^'>° J -L JR. 

Likewise, we want to establish conditions for when 

drlf f{C,;r)dc) =([ drf{C,;r)dc) , 

XJM. I y XJR I y 

which is equivalent to showing 

lim lim ^ j [ \ fr,h{C, y, r) - fr{x, y, t)| d^dy = 0, 

where 

We assume that f{^,y,T), drf{C,y,T) G L^{M. x [—L,L]) on some interval 
r G (ti,T2). Then we have, using dominated convergence, for any L > 0, 
that ^ 

l^^r.^ / IfrM^^y^T) - fr{x,y,T)\dCdy = 0. 



Define the function f{L, h) such that 



fiL,h) = ^ [ [ \fT,h{^,y,T)- fr{x,y,T)\ 

J-L JR 



didy. 
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We assume supi-^i^* f{L, h) < oo for some L* > 0. Then, for all e > 0, there 
must be some value L such that 

limsup-^- / / \fr,h{C,y,T) - fr{x,y,T)\d(,dy < f{L,h) + e. 

L-5-00 ^-l^ J-L Jr 

Taking the limit on both sides as /i — )■ 0, and noting that e > is arbitrary 
shows 

linilimsup-^ / / \fr,h{^,y,T)-fr{x,y,T)\d^dy = 0. 

L-s>oo '^^ J-L JR 

Since we are only working with positive quantities, this shows that 

lim lim / / |/r,h(C, y, - /r(a;, y, r)| d^dy = 0, 

and the result is established. Thus we have found a class of functions for 
which the analytic operations performed throughout the paper are valid. 
We believe the web solutions and their perturbations satisfy the conditions 
listed above, but it is beyond the scope of the paper to prove as such. 



Conservation of Energy for KP Web Solutions 



In this section, we show the condition (17) holds for web solutions. We also 



establish a number of properties about the web-solutions with regards to the 



average used in defining the energy. Again we point out that if (17) holds, 
then from 



WWrdS^ 



wK{w)d^ 



I Wdt ^Wyyd^ 

Jr 



where we have assumed perfect derivatives in ^ cancel, dr£{w) = 0. To get 



the derivative of the energy to vanish then, we must prove (17) for u = w. 
To do this, using integration by parts, we find 

fL 



1 



wdt ^Wyydyd^ = — wdc ^ 



1 



L 



y=L 

-L 



di 



(33) 



yd^ ^Wydyd^. 



Since we have that Wyd^ ^Wy = ^d^{d^ ^'^yf' ^'^^ expect that {d^ ^ 

y=L 



0, cf. [7], and 4- wdc 



as |L| — 7- oo, then the result should follow. 



To prove (17) rigorously, we begin by showing that w is always positive. 



since this allows us to introduce notation and simplify some computations. 
To show w > 0, note that w can be written as 



w 
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and that we can write in the form Q = det{BEK), where B is a N x M 
matrix with entries bij, E is an M x M diagonal matrix with Ejj = e^^ , 

and i^T is an M X matrix with entries Kij = fcp {cf. [llj). Using the 
Cauchy-Binet Theorem, 17 can be written as 



n= Y det{B<N>,s)det{{EK) 



S,<N>), 



where (M) = {1, 2, • • • , M} and ^^^^^ represents all possible permuta- 
tions of the M numbers in (M). The symbol i?<Ar>,5 denotes the N x N 
matrix formed from the S = {S'(l), 5(2), • • • ,S{N)} columns of B. Like- 
wise, {EK)s,<:N> is the N x N matrix formed from the S rows of the matrix 
EK. By construction it is assumed that det(i?<Ar>,5) > 0, and that not all 
of these terms can be identically zero. Likewise, for a given S, it follows 
that 



det{{EK)s,<N>) = exp ("^ Osn)] H 

\l=l J l<i<j<_ 



iks(j) - ks{i)). 

<j<N 

Therefore J7 > since ks(j) > ks(^i-^ when j > i. Letting 

As = det{B<N>,s) n ^^sij) - ks(i)), 

l<i<j<N 

it is straightforward to show that Q^^fi — (fi^)^ is equal to 

(N \ N N 

j=l J 1=1 n=l 

In the sums, if we have the pair {S,S'), then we must also have the pair 
(5', S), and combining these terms gives 

As As' exp 1^ 9s(j) + 9s'ij)^ {ks{i) - , 

and therefore w > 0. 

Since Q is smooth in all its arguments, w is as well, and we then have 



rL2 

/ w{ty,T,T)d^ = 2d^lnn 

Jli 



We need to determine the behavior of rig /fi as L2 — >■ oo and Li — )• —00. In 
the case that L2 — )• 00, there is an S, say Smax = {ku-N+i, • • • , kpi}, such 
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that for C > 1, 

' N \ IN 



exp I Y,^s^a.{3) > exp Y.^s(j) 



We write Vl^/n as 

Z=l 

where we have divided through by A^^^^ exp {Y1^=i ^Smax(j)) *o isolate the 
leading order term. An identical argument shows that 

where Smin = {1; 2, • • • , A^}. Therefore, we have shown that 

\i=i 1=1 J 

This result also establishes that the average (/jj wd^)^ is well defined. Using 
the inequality 

/ / w'^{^,y,T,T)dCdy <sup\w\ ^ [ [ w{S,,y,T,T)dS,dy 
yields, noting that supig2 \w\ < oo, 

limsup— — / / w'^{£,, y, T, T)d^dy < oo. 

L-5-00 2L JiR 

However, establishing that the limit exists, or that £{w) is rigorously well 
defined for any web solution is technically more demanding and beyond the 
scope of this paper. Therefore, we consider the energy £{w) formally. 

Returning to evaluating (33), for the second term on the right hand side, 
using Wyd^^Wy = \d^{d'^^Wyf, we see we need to evaluate 



lim [dj^w,, 



i=-Li 



We have that 



dc ^Wy = 2dy^ In r2 — lim dy^ In — lim dy^ In fi. 
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Given dy^ In Vt 



we write Vt^yQ. — VLfVLy as 



N 



N 



(i) 



S,S' 



Thus, in the numerator of , terms in the sum where S = S 

cancel, which does not happen in the denominator. Thus, as ^ — >• c«, the 
dominant behavior of the denominator, , , which is exponentiaUy 

larger than any term in the numerator. A similar argument holds for ^ — )• 
— oo, and thus lim|g|_^oo dy^\nVL = 0, and therefore, for any L > 0, 



Wyd^ ^Wydyd^ = 0. 



Finally, we now need to show that 



lim 

L->oo 



- Wd^ Wy 



y=L 



di = 0. 



It is straightforward to show that 



^Wy 



is uniformly bounded above 



throughout the i — y plane, say 
Thus 



^Wy 



< C, so that 



< Cw. 



- Wd^ Wy 



y=L 



di 



< y(9^1nl7(^,L,r,r) + aglnO(^,-L,r,r)) 



which immediately gives the estimate 



- Wd^ Wy 



y=L 

-L 



2C 



N 



N 



1=1 



and therefore we obtain the desired result (17). 

Calculation of the pseudo-spectral representation of d^^ 

We take the number of interpolation points in the pseudo-spectral method 
to be Nt, which is assumed to be an even number, i.e. Nj- = 2N. Then 
each interpolation point is given by 



2L. 
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with m = 0, • • • , Nt — 1 and n = 0, • • • , — 1. We define the coefficients 
aji{T) = (— l)''^-'aj7(r), so that at the point {^m,yn), using (30), we have 

LIT ^ ^ ^ ' 



ivr / 



2-KijnjNx 

^ I 

j=~N+l /=-Ar+l,/^0 

TV 

j=-N+l 

We then, for integers j and /, compute the discrete inverse Fourier transform 

Nt-1 Nt-1 



which gives 



^ m=0 n=0 



I "jOV ^ t -2-Kiml/NT 

m=0 



If / = 0, the last sum in this expression becomes —L^. Otherwise, if / 7^ 0, 
we have 

Nt-1 j Nt-1 

^-27Timl/NT ^ -^d- ^-2mml/NT 

m=0 m=0 

Le^ ( - 1 \ 



ITT ' I g-2TTil/NT _ 1 y 

where we treat / as a continuous parameter. After differentiating, we get 
the final expression 



i=-JV+l,/7^0 



which is (31) after relabeling indices. 
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